function hfile2 = vmr_CleanVMR(hfile, opts)
% VMR::CleanVMR  - try automatic cleaning (experimental !)
%
% FORMAT:       cleaned = vmr.CleanVMR([opts])
%
% Input fields:
%
%       opts        1x1 struct with optional settings
%        .framed    if true, assumes the VMR is already framed
%
% Output fields:
%
%       cleaned     new, cleaned VMR object

% Version:  v0.7g
% Build:    9100720
% Date:     Oct-07 2009, 8:48 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmr')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if nargin < 2 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'framed') || ...
   ~islogical(opts.framed) || ...
    isempty(opts.framed)
    opts.framed = false;
else
    opts.framed = opts.framed(1);
end

% copy vmr and get double precision data
hfile2 = aft_CopyObject(hfile);
sfile2 = bvqxfile_getscont(hfile2.L);
sfile2.F = '';
if ~bc.VMR8bit || ...
    isempty(bc.VMRData16) || ...
    numel(size(bc.VMRData16)) ~= numel(size(bc.VMRData)) || ...
    any(size(bc.VMRData16) ~= size(bc.VMRData))
    cvmr = class(bc.VMRData);
    vmrd = double(bc.VMRData(:, :, :));
else
    cvmr = 'uint16';
    vmrd = double(bc.VMRData16(:, :, :));
end
isz = size(vmrd);

% if not framed, inspect exterior planes first
if ~opts.framed
    pxy1 = vmrd(:, :, 1);
    pxye = vmrd(:, :, end);
    pxz1 = squeeze(vmrd(:, 1, :));
    pxze = squeeze(vmrd(:, end, :));
    pyz1 = squeeze(vmrd(1, :, :));
    pyze = squeeze(vmrd(end, :, :));
    
    % build mean / std of each plane
    mxy1 = mean(pxy1(:));
    mxye = mean(pxye(:));
    mxz1 = mean(pxz1(:));
    mxze = mean(pxze(:));
    myz1 = mean(pyz1(:));
    myze = mean(pyze(:));
    sxy1 = std(pxy1(:));
    sxye = std(pxye(:));
    sxz1 = std(pxz1(:));
    sxze = std(pxze(:));
    syz1 = std(pyz1(:));
    syze = std(pyze(:));
    
    % find planes for which the mean is around the mean of means
    dmea = mean([mxy1, mxye, mxz1, mxze, myz1, myze]);
    mstd = 0.5 + mean([sxy1, sxye, sxz1, sxze, syz1, syze]) / 1.5;
    tv = [];
    if abs(mxy1 - dmea) < mstd
        tv = [tv; pxy1(:)];
    end
    if abs(mxye - dmea) < mstd
        tv = [tv; pxye(:)];
    end
    if abs(mxz1 - dmea) < mstd
        tv = [tv; pxz1(:)];
    end
    if abs(mxze - dmea) < mstd
        tv = [tv; pxze(:)];
    end
    if abs(myz1 - dmea) < mstd
        tv = [tv; pyz1(:)];
    end
    if abs(myze - dmea) < mstd
        tv = [tv; pyze(:)];
    end
    
    % any background found
    if ~isempty(tv)
        
        % set all voxels up to mean + 0.5 * std + 1 auf 0
        coff = mean(tv) + 0.5 * std(tv) + 1;
        vmrd(vmrd <= coff) = 0;
    else
        coff = 0.5;
    end
    
    % find remaining data
    nxc = [isz(2) * isz(3), 1];
    for x1 = 1:isz(1)
        if any(reshape(vmrd(x1, :, :), nxc) > 0)
            break;
        end
    end
    for xe = isz(1):-1:x1
        if any(reshape(vmrd(xe, :, :), nxc) > 0)
            break;
        end
    end
    sxc = x1:xe;
    nyc = [numel(sxc) * isz(3), 1];
    for y1 = 1:isz(2)
        if any(reshape(vmrd(sxc, y1, :), nyc) > 0)
            break;
        end
    end
    for ye = isz(2):-1:y1
        if any(reshape(vmrd(sxc, ye, :), nyc) > 0)
            break;
        end
    end
    syc = y1:ye;
    nzc = [numel(sxc) * numel(syc), 1];
    for z1 = 1:isz(3)
        if any(reshape(vmrd(sxc, syc, z1), nzc) > 0)
            break;
        end
    end
    for ze = isz(3):-1:z1
        if any(reshape(vmrd(sxc, syc, ze), nzc) > 0)
            break;
        end
    end
    ioff = [x1, y1, z1];
    isz = [xe, ye, ze];
    vmrd = vmrd(x1:xe, y1:ye, z1:ze);
else
    
    % set offsets to 1
    ioff = [1, 1, 1];
    coff = 0.5;
end

% smooth remainder
svmrd = smoothdata3(vmrd, [3, 3, 3]);

% build histogram
[hn, hx] = hist(svmrd(svmrd(:) > coff), 100);
ha = sum(hn);

% find peak in histogram (supposedly with remainder of noise)
[mxn, mxi] = max(hn);

% check surroundings
if mxi > 1 && ...
    mxi < 50 && ...
    hn(mxi - 1) < (0.75 * mxn) && ...
    hn(mxi + 1) < (0.75 * mxn) && ...
    hn(mxi + 2) < (0.33 * mxn)

    % find threshold adaptively
    tx = mxn;
    ha = (ha - (hn(mxi - 1) + mxn + hn(mxi + 1))) / 97;
    for tc = mxi:100
        if hn(tc) > tx || ...
            hn(tc) < ha
            break;
        end
        tx = hn(tc);
    end
    tx = hx(tc) + 3;
    
    % remove upto threshold + 3
    vmrd(vmrd < tx) = 0;
    tx = tx / 3;
    
    % smooth again and repeat once more to remove high intesity speckles
    clear svmrd;
    svmrd = smoothdata3(vmrd, [3, 3, 3]);
    vmrd(svmrd < tx) = 0;
    tx = tx / 3;
    clear svmrd;
    svmrd = smoothdata3(vmrd, [3, 3, 3]);
    coff = tx;
end

% clear vmrd
vmrd(svmrd <= coff) = 0;

% if biggest cluster from center, only keep that
vmrs = (vmrd > 0);
numv = sum(vmrs(:));
ccr = ceil(isz / 2);
vmrs = floodfill3(vmrs, ccr(1), ccr(2), ccr(3));
if sum(vmrs(:)) > (numv / 2)
    vmrd(~vmrs) = 0;
end

% back to old class
if strcmpi(cvmr, 'uint8')
    vmrd = uint8(vmrd);
else
    vmrd = uint16(vmrd);
end

% resolve transio first
if istransio(sfile2.C.VMRData)
    sfile2.C.VMRData = ...
        reshape(bc.VMRData(:), size(bc.VMRData));
end
if istransio(sfile2.C.VMRData16)
    sfile2.C.VMRData16 = ...
        reshape(bc.VMRData16(:), size(bc.VMRData16));
end

% where to put this
if ~sfile2.C.VMR8bit || ...
    isa(vmrd, 'uint8')
    sfile2.C.VMRData(:) = 0;
    sfile2.C.VMRData(ioff(1):isz(1), ioff(2):isz(2), ioff(3):isz(3)) = vmrd;
else
    sfile2.C.VMRData16(:) = 0;
    sfile2.C.VMRData16(ioff(1):isz(1), ioff(2):isz(2), ioff(3):isz(3)) = vmrd;
    sfile2.C.VMRData(sfile2.C.VMRData16 == 0) = 0;
end

% put content back
bvqxfile_setscont(hfile2.L, sfile2);
